// PowerSpectrum.cpp: implementation of the PowerSpectrum class.
//
//////////////////////////////////////////////////////////////////////

//#include "stdafx.h"
#include<stdio.h>
#include<string.h>
#include "PowerSpectrum.h"
#include<iostream.h>
//:wq
//#include<math.h>
#include<complex>
using namespace :: std;

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

PowerSpectrum::PowerSpectrum()
{
	s=NULL;
	x=NULL;
	strCRegion=NULL;

}

PowerSpectrum::~PowerSpectrum()
{
	delete []s;
	delete []x;
    delete []strCRegion;

}

void PowerSpectrum::setSequence(const char * seq)
{

delete []strCRegion;

strCRegion=new char[ strlen(seq)+1 ];

memset(strCRegion, '\0' , strlen(seq)+1 ) ;

strcpy(strCRegion,seq);

}


int PowerSpectrum::doPowerSpectrumAnalysis() 
{


const double PI = 3.1415926535897932384626433832795;
//public static final double PI ;


complex<double>  i1(0.0,1.0);
//complex<double>  i1 = sqrt(-1);

long lenStr = strlen(strCRegion);

delete []s;
delete []x;


cout<<"Coding Region Length:"<<lenStr<<endl;


s	= new double[lenStr/2 + 1];
x	= new double[lenStr/2 + 1];

//temporary check; 
x1  =  new double[lenStr/2 + 1];
xp  =  new double[lenStr/2 + 1];
xp1  =  new double[lenStr/2 + 1];
xr  =  new double[lenStr/2 + 1];
xr1  =  new double[lenStr/2 + 1];
 
double *sf  =  new  double[lenStr/2 + 1]; 

char nuclt[5]= "ATGC";
int const lenWindow = (lenStr/3)*3;



double winsq		= (1.0/lenWindow)+(1.0/(lenWindow*lenWindow));
long  wincube	    = lenWindow*lenWindow*lenWindow;


int countA =0;
int countC =0;
int countG =0;
int countT= 0;

complex<double> s1(0.0,0.0);
double ss2=0.0;
double sq=0.0;
double s2=0.0;
 //temporary check
double p1=0.0;
double p2=0.0;
double p3=0.0;

double sbar =0.0;
double p=0.0;

	
//ofstream out("power");

	for(long i =0; i < lenWindow; i++ )
	{

	if( strCRegion[i] == 'A') 
		{ countA++; continue;}
  	if( strCRegion[i] == 'G') 
		{ countG++; continue;}
  	if( strCRegion[i] == 'C') 
		{ countC++; continue;}
  	if( strCRegion[i] == 'T') 
		{ countT++; continue;}

	} 

//out<<"Count A  :"<<countA <<"    countT :"<<countT <<"countG   :"<<countG <<"countC   :"<<countC <<endl;

 
sq = 0.0+ pow(countA,2)+pow(countG,2) + pow(countT,2) + pow(countC,2);

double f =0.0;

for(long k = 1 ; k <=lenWindow/2  ; k++ )
{

	s2		=0.0;
	sbar	=0.0;
	p		=0.0;
	f= ((1.0)*k)/lenWindow;

// for each Nucleotide. A T G C


	for(int i = 0; i <4 ; i ++ )
	{
		s1= complex<double>(0.0,0.0);


		for (int bi=0; bi<lenWindow; bi++)
		{
				if(strCRegion[bi] == nuclt[i] )
					s1= s1+exp(i1*2.0*PI*f*(bi+1.0));
		}

		s2 = s2+ (abs(s1)*abs(s1));
		//ss2=ss2  + (1/lenWindow*lenWindow)*(abs(s1)*abs(s1));

	}

	//store the value of sf for each freq.
	sf[k-1]=s2;

//	to calculater power specturm and code region 

		sq   = sq/wincube;
		sbar =(winsq-sq);
		sq   = sbar * (lenWindow/2.0);

// to calcuater signal to noise ratio

		s2 = s2/(sq*lenWindow*lenWindow);

		p= s2/sbar;
  //temporary check;
        //p1=(p*lenWindow*sbar)/2;

	
		//p2=p1*2;
		p2=p/(lenWindow*sbar);
		p3=2*p2;

	s[k-1]=f;
	x[k-1]=p;
	xp[k-1]=p2;
	xp1[k-1]=p3;


	//temporary check
	//x1[k-1]=p1;
	//x2[k-1]=p2;
	//x3[k-1]=p3;


//	out<<"\n Frequency:"<<f <<" Power  : "<<p;

} // outer for loop

//calculating the sbar value

 double sbar1=0.00;
 double sbar2=0.00;

  for(k = 1 ; k <=lenWindow/2  ; k++ )
  {
	  sbar1=sbar1+sf[k-1];

  }
  sbar1=(sbar1*2)/lenWindow;
  sbar2=sbar1/(lenWindow*lenWindow);
  //calculating the sf/sbar1
  

  for(k = 1 ; k <=lenWindow/2  ; k++ )
  {
	   x1[k-1]=sf[k-1]/sbar1;

	  xr[k-1]=x1[k-1]/(lenWindow*sbar2);
	  
	  xr1[k-1]=xr[k-1]*2;

  }



int numPeaks=0;
int genePeak=0;

for( k = 1 ; k <= (lenWindow/2)  ; k++ )
{

	if( x[k-1] > 4.0)
	{
	   //if(  s[k-1] == 1.0/3.0  )] //added by Sreenu on mandars suggestion on 19/04/2002
	 if( ( s[k-1] >= 0.32 ) && ( s[k-1] <= 0.34 ))
             genePeak=1;
	   else
		   numPeaks++;
	}

}


//delete []s	  ;
//delete []x	  ;
//s =NULL  ;
//x =NULL  ;


if(genePeak==1)
{  
	if(numPeaks==0)
		 return 1;
	 else
		 return 2;
}
else
{
	return 0;
}



}  // end of member function.



/*PowerSpectrum& PowerSpectrum::operator =(const PowerSpectrum &ps)
{

}*/
